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ABSTRACT 


Standard theoretical methods of analysis which work well for seem- 
ingly more complex problems fail to predict the experimentally observed 
instability of fully developed, incompressible pipe flow at any Reynolds 
number. Past research by Harrison and Arnold on the stability of pipe 
flow yielded erroneous results due to errors in the setup of the problem 


and formulation of the boundary conditions at the axis. 


A revised theory with particular attention to the rather complex 
boundary conditions at the axis has recently been developed. Improved 
finite differencing techniques with consistent fourth order truncation 
error were also used to approximate the governing differential equations. 


Numerical results for angular wave numbers zero and six show that 


the flow is stable at all Reynolds numbers. Results for angular wave 


number one contain instabilities at all Reynolds numbers for small 


values of the axial wave number. These results are tabulated, plotted, 


and discussed in detail in this mK 
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I. INTRODUCTION 


Osborne Reynolds [Ref. 1] conducted his classical experiments on the 
transition from laminar to turbulent flow of fluids in circular pipes 
nearly 100 years ago. Based on his experiments the critical Reynolds 
number for pipe flow was found to be 2300. An analytical solution for 
the problem of predicting instabilities of fully developed, three dimen- 
sional, incompressible flow of constant viscosity fluids in pipes has 
been pursued actively since then. Several approaches have been taken in 
solving the inherently nonlinear Navier-Stokes equations, which along 
with the continuity equation, govern the behavior of fluid flow in pipes 
of circular cross-section. 

Salwen and Grosch (Ref. 2] studied pipe flow with purely sinusoidal 
streamwise (axial) perturbations. Infinitesimal velocity and pressure 
perturbations, which were explicit functions of time and not complex or 
exponential in form, were introduced into the flow field. Numerical 
calculations were carried out at various angular wave numbers and it was 
concluded that the flow was stable for all axial wave numbers and 
Reynolds numbers. Perturbations with exponential growth in space but a 
purely sinusoidal variation in time were explored by Garg and Rouleau 
(Ref. 3] and found to be stable. Gill (Ref. 4] studied combinations of 
exponential growth in space and in time using power series analysis and 


concluded that these flows were al] stable. 
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Because of the general agreement that pipe flow is stable to infini- 
tesimal disturbances, two other approaches to the problem of predicting 
a instabilities at the experimentally observed critical Reynolds number 
i have been pursued. First, investigations by Davy and Drazin (Ref. 5] 
confirmed that pipe Poiseuille flow was stable at all Reynolds numbers 
for infinitesinal disturbances that were both temporally and axially 
complex and exponential in nature. However, despite the fact that the 
governing equations must necessarily remain nonlinear with the introduc- 
tion of disturbances of finite amplitude, they concluded that the flow 
is unstable with respect to finite disturbances. Similar theoretical 
examinations of plane Poiseuille flow between two parallel plates led 
McIntire and Lin [Ref. 6] to this same conclusion. Second, it was 
postulated by Huang and Chen [Ref. 7] and Leite (Ref. 8] that the origin 
of flow instabilities occur in the entrance region of the pipe, where 
the flow has not yet become fully developed. Both theoretical and 
experimental studies have been done to support these hypotheses. Garg 
(Ref. 9] also found that flow instabilities existed near the entry 
region for developing pipe flow. He concluded that results using the 
Hornbeck velocity profile more nearly approximated experimentally deter- 
mined values of the critical Reynolds number. While these investi- 
gations have shown instabilities to exist in pipe flow, a completely 
generalized solution to the linearized problem of fully developed pipe 
flow has never been accomplished. 

Recently, a more general theory consisting of perturbations which 


have a fully complex exponential form with respect to time and the axial 
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coordinate, a purely imaginary exponential form in the angular coordi- 

i nate and are three dimensional in nature have been explored. Develop- 

. ment of this theory by Harrison (Ref. 10] and further numerical in- 

vestigations by Arnold (Ref. 11] failed to produce conclusive results 

that instabilities exist in fully developed pipe flow. Errors in the 

setup of the problem and inadequate formulation of the boundary condi- 

tions at the pipe axis contributed to these erroneous results. Arnold, 

N however, was on the right track in his research by incorporating newly 

formulated boundary conditions at the axis developed by Gawain (Ref. 12]. 

Arnold's work, for the most part, remains valid except that he 

overlooked a small but significant detail that resulted in unclear 

: definitions of stability and instability. He used a complex axial wave 

: number a, of the form a = dp + iay, whereas the correct version should 

; simply be a purely imaginary quantity of the form ia. Thus, Arnold 

introduced an extra degree of freedom, namely the fictitious and in- 

: correct quantity Gp: This accounts for the erroneous instabilities 
which he obtained in his numerical investigation of pipe flow. 

Further work by Gawain (Ref. 13] refined the boundary conditions at 
the axis and incorporated the previously defined imaginary axial wave 
number, in exponential form, which represents only sinusoidal oscilla- 
tions with respect to the axial coordinate. Using these advancements in 
the linearized theory, Gawain showed an asymptotic trend toward neutral 
stability of the flow with increasing Reynolds numbers for angular wave 
numer, n = 0. This paper duplicates and confirms the results for 
n=0, that the flow is stable at ail Reynolds numbers and axial wave 


numbers. It is further shown that fully developed pipe flow for n = 6 
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is also stable and shows the same trends as for n= 0. Instabilities 
are shown to exist in the numerical analysis of the vorticity transport 
equations for n = 1. The paradox that now appears is that for angular 
wave numbers n = 0 and 6 the flow is apparently stable at all Reynolds 
numbers while for n = 1 the flow is apparently unstable at all Reynolds 
numbers. Neither of these theoretical results is consistent with the 
known experimental fact that there exists a critical Reynolds number 
below which pipe flow is stable and above which it is unstable. However, 
the new result for n = 1, while it can hardly be called correct, is at 
least encouraging in that it shows for the first time that instabilities 
can in fact exist in the solution of the linearized vorticity transport 
equations for fully developed pipe flow. 

In addition to implementing the improved axis boundary conditions 
and the advancement in the linearized theory, improved finite differ- 
encing techniques were used to approximate the vorticity transport 
equations along a uniform radial mesh. Highly accurate numerical solu- 
tions for the previously mentioned angular wave numbers, n = 0, l, 
and 6, were obtained using double precision, complex numbers. Addition- 
ally, finite difference equations with consistent fourth order trunc- 
ation error were used throughout. Details of the development of the 
vorticity transport equations and the central and non-central finite 
difference equations are discussed jin later chapters. Systematic and 
extensive calculations remain to be accomplished for angular wave 
numbers, n = 2, 3, 4, 5, and n > 6. The results of these numerical 


analyses should prove to be fruitful in establishing a theory which 


adequately explains the experimentally observed fact of the onset of 
flow instabilities at the critical Reynolds number. 
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II. THE VORTICITY TRANSPORT EQUATION ~° 


A. BASIC THEORY 


The governing equations for laminar flow of fluids in a circular 


pipe are derived from Newton's second law of motion, 


F=uma (2-1) 


a and the conservation of mass principle through an infinitesimal contro] 


volume. 


Equation (2-1) results in the familiar Navier-Stokes equation, 


which is a vector equation in three cylindrical coordinates and in time. 
The equation for conservation of mass results in the continuity equation 


for incompressible flow and is given by, 


v-V=0 (2-2) 


For the case at hand, the following assumptions are made, 


1. 


Ress ara, 


“2 7 Pp wD 


The fluid is incompressible, p = constant 

The fluid has constant viscosity, p = constant 
The flow is fully developed, U = f(r) only 

The flow is three dimensional in nature 

The mean flow is steady, av/at = 0 

The perturbation flow is unsteady, av/at # 0 


The effects of body forces are negligible 


The resulting Navier-Stokes equation can be stated as the sum of the 


forces per unit mass is equal to the acceleration or, 


Bat ie es Rage Oe Ec ae 
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z mane acceleration (2~3) 


The forces acting on the fluid are pressure forces and viscous forces. 
Body forces and hydrostatic forces balance out and are not a factor 
here. 

Equation (2-3) in its dimensional form can be expressed as, 


- ve Eis (ie) vee te a | (2-4) 


where the starred quantities are fully dimensional. 
It should be noted here that the equations presented through equa- 
tion (2-18) represent the mean flow. Equation (2-4) can also be 


expressed in a nondimensional form as, 


-vp+lL gyal (2-5) 


where Re is the Reynolds number and is defined as, 


a 
(=f ] 


Re 


u 
<| 


(2-6) 


This definition of Reynolds number is based on the pipe radius R as the 
characteristic length, the mean volumetric velocity U as the character- 
istic velocity, and the kinematic viscosity v which is constant. The 
three quantities on the right side of equation (2-6) are dimensional, 
but Re is a dimensionless quantity. On this basis, the critical 
Reynolds number for transition from laminar to turbulent flow becomes 


1150 (vice the value of 2300 which is obtained if the Reynolds number is 
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defined in teres of pipe diameter). All subsequent equations presented 
in this paper are in nondimensional fora. 

Returning to the derivation of the mean flow, the acceleration term 
or substantial derivative, dV/dt can be expanded to the form, 


Hog. oy. 2 (2-7) 
but (V - 7)¥ can be further expanded by a well known vector identity to, 
@- wis (F) - Fx xd) (2-8) 


Additionally, V2V from equation (2-5) can be expanded to, 
vzV¥=ev(9-V)-7x (Vx ¥) (2-9) 


but V - V = 0, since continuity must be satisfied and equaton (2-9) 


becomes, 
wz¥=-v0x (Vx ¥) (2-10) 


With appropriate substitutions, the final form of the Navier-Stokes 


equation becomes, 
7 rae | ay U2 _% 2 av = 
VP we VX OxV a0) Vx (Vx V) + 5 (2-11) 
A more convenient form of equation (2-11) is, 


-V (P+ (Py) shox (vx - Vx xt) +8 (2-12) 
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As previously discussed, al] the quantities in the equations have 
been nondimensionalized. Dimensionless cylindrical coordinates x, r, 
= and © are used with Can é. and % denoting the unit vectors in the three 
cS coordinate directions, respectively. It is well known that the velocity 
_ profile for fully developed, three dimensional, incompressible flow in 

- circular pipes can be expressed analytically by, 
a U=2(1- r2) (2-13) 
R} This function has the shape of a paraboloid of revolution as shown in 
3 Figure 2-1. 
i 
i 
¥ 


Figure 2-1. Velocity Profile of Fluid Flow in a Pipe 


The nondimensional mean velocity vector V(r) of the flow is given by, 


Wor) =U é, =2(1- r2) é (2-14) 
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A second vector quantity a, the mean vorticity, is defined as the curl 
of the mean velocity vector and is given by, 


GavxVeara, (2-15) 
which will be used later. 
F B. THE VORTICTY TRANSPORT EQUATION 


The vorticity transport equation is now obtained by taking the cur! 
of the Navier-Stokes equation (2-12), resulting in, 


: -vxv P+ Oy) shy x ex (vx) - 9% Mx (Wx MH) + XD 


(2-16) 


The primary reason for taking the curl of the Navier-Stokes equation 
is to eliminate the unknown scalar on the left side of equation (2-16). 
Since the curl of the gradient of a scalar is equal to zero, equation 


- (2-6) reduces to, 
ae VX (VX (VX 9D) - 0x x (Vx H) + AREY 2g (2-17) 


Equation (2-17) can be simplified by substituting the mean vorticity 


vector relation equation (2-15) into equation (2-17) resulting in, 


at 


he VX (Wx dy -0x Wx %) +B 


= 0 (2-18) 
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In order to analyze this flow field and in particular to determine 
the transition from laminar to turbulent flow, a disturbance of sual] 
amplitude is introduced into the flow. It is assumed that the resulting 
flow is made up of the steady state, laminar flow and a small perturba- 
tion flow superimposed on the laminar flow. The total velocity and 
vorticity respectively become, 


2 Ge8 (2-19) 
and 
4°23 +0 (2-20) 


where v and w are the velocity perturbations and vorticity perturbations, 


respectively. It is of course a necessary requirement that the contin- 


uity equation for the total flow be satisfied here as well, 
v-Ws0-¥ + 9+ V=0 (2-21) 


In view of this relation and of equation (2-2), the perturbation flow 


independently satisfies the continuity condition and 


OT By 


The introduction of the perturbation quantities into the vorticity 


y 


transport equation is accomplished by substituting V’ for V and &° for 4 


into equation (2-18). This results in, 


ee oD ede ala 


Lox (vx +a) - 0x G+) x Gedy) + Geo 24 


(2-23) 


3 24 


~ we Mae 
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which is the equation for the total or resultant flow field. Upon ex- 
panding equation (2-23), applying the mean steady flow assumption, 
av/at = 0, and subtracting equation (2-18) from it, equation (2-23) 


becomes, 


ol®s 


ts Vx (Vx) - 7x (xe) - xv) + (x )] + Me 0 


(2-24) 


Now equation (2-23) is linear in the perturbation quantities except for 
the second order term, (Vv x w). Since the perturbations are assumed to 


be small, this term can be neglected. 


The perturbation vorticity is defined as, 
w=EUXx (2-25) 


which is analogous to equation (2-15) for the mean flow. Making this 
substitution in equation (2-24) yields the linearized vorticity trans- 
port equation, 

1 av 


fe VX (VX (VX V)) - 0X We (x V)) FOX (VX Hx Y)) + OK He0 


(2-26) 


As previously mentioned, the continuity equation (2-22) for the per- 
turbation velocity must be satisfied. A convenient way to insure this 


is to express v in terms of a velocity vector potential function W as 


[ foi lows, 


v=uxwW (2-27) 


It can be verified that if v is defined in this way, equation (2-22) is 
satisfied identically for any arbitrary vector function W. 

Instead of proceeding directly with the real vector functions Vv and 
W, it is advantageous to work with the Fourier transforms of Vv and W, 
which are in general complex functions. This is the form which is used 
here. It can be shown, because of the linearity of the basic equation, 
that if a solution can be found for the Fourier transforms, then the 
corresponding real vector functions will satisfy the basic equation as 
well. A detailed explanation of these functions was presented more 


elegantly by Gawain [Ref. 13]. The vector potential function W, now 


representing the Fourier transform, can be expressed as, 
W= (F(r) @, + G(r) @. + H(r) é] e* (2-28) 


where X is a convenient abbreviation for the complex quantity, 
X = jax + ind + yt (2-29) 


The axial wave number a is a real quantity greater than or equal to zero 
which in the exponential form represents purely sinusoidal axial vari- 
ations of v. The angular wave number n is a positive integer and repre- 
sents sinusoidal angular variations of Vv which are periodic with respect 
to coordinate 6. The complex frequency y contained in the exponent is 


of the form, 
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okt = g(%y * ¥pt = Q¥gt (cos yp t+ isin y, t) (2-30) 


The quantity Yep clearly represents the exponential time rate of growth 
or decay of the perturbation amplitude. The basic criteria for hydro- 
dynamic stability can be defined in terms of Yp- Positive values of Yp 
represent growth and signify flow instability while negative values of 
Yp represent decay and signify flow stability. 

The vorticity transport equation can now be expressed in terms of 
the velocity vector, V and the complex vector potential, W. By making 
the substitution for v, equation (2-26) now becomes, 


pi 


te vx (Vx (Vx (Vx W)) - 0x (Wx (9 x (0 x WD) 


tox (Vx Hx (Vx w)) + 0% (Vx BM <0 (2-31) 


In convenient shorthand notation, equation (2-31) becomes, 


Equation (2-31) = f(r)e* = [r.(r) é. + T(r) e, *ra(r) @je~=0 


(2-32) 
and can be expressed in matrix notation as, 
Cy 0 
rs te =4 0? (2-33) 
ley S 
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Therefore, the three components of the vorticity transport equation are 


equivalent to three simultaneous scalar equations, namely, 


r. 20 (2-34) 


Of these three equations, only two are actually independent since, 
as can be seen from equation (2-31), f e*% is itself the curl of another 
vector. Therefore f e* must necessarily be nondivergent. It can be 


shown that the components of fr satisfy the identity, V - ‘(a e X) = 0, or, 


iol, +loory + 2) T, =0 (2-35) 


where D is the symbol for the differential operator d/dr. Two linearly 
independent equations can now be formulated from equations (2-34) given 
the constraint of equation (2-35). This can be reduced to some appro- 
priate linear combination of two of the equations which is set equal to 
zero and to the third equation which is also set equal to zero. This 
results in a system of two coupled, fourth order differential equations 
in three unknowns. Recal] that the vector potential function W is 
expressed in terms of F(r), G(r), and H(r). It can be shown that one of 
these variables is redundant and may be eliminated without any loss of 
generality in the solution of the equations. For the present case, the 


solution assumes its most convenient form if, 


F(r) = 0 (2-36) 


This convention has therefore been adopted and now G(r) and H(r) become 
the two coupled eigenfunctions of the governing fourth order differen- 
tial equations. Returning to the problem of selecting the proper linear 
combination of equations comprising f, it is important to choose this 
combination in an optimum fashion by analyzing the algebraic structure 
of these equations. Table 2-1 is taken from reference [13], corrected, 
and duplicated here for that purpose. Notice that equations (1) and (3) 
in the table are third order in G. If the particular linear combination 
of these equations is chosen as appears in equation (4), the result 
reduces to an equation of second order in G. This does not affect the 


order of the result since the function H remains fourth order. 


Table 2-1. Properties of the Vorticity Transport Equations 


HIGHEST 
TERM OF TERM OF NEGATIVE 
HIGHEST HIGHEST POWER 
EQUATION ORDER IN G ORDER IN H OF r 


ia ps 
Re 0°G 


- str (r)+rg(r)=0 
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The forsulation of the problem is simplified greatly by choosing two 
independent vorticity transport equations in the fora, 


ra 2-0 (2°37) 
and 
= = ie 
he (2-38) 


C. THE VORTICITY TRANSPORT MATRIX EQUATION 


While the fourth order linear vector operations and the associated 
algebra becomes very lengthy and rather tedious, the investigator at- 
tempting to expand the final form of the linearized vorticity transport 
equation (2-31) should pay sufficient attention to these details. 
Recall that the curl operation on a vector in cylindrical coordinates 


takes the form, 


1? 1+ > 
re’ 7 °F &9 
-/2 @& 2 : 
a ean me or 36 (2-39) 
0 Ge* rHe* 


After all the nabla (V) operations are complete, the final equation wil] 
be in the form of a vector equation as defined in equation (2-32). The 
first vorticity transport equation consists of the coefficient repre- 


sented by rh and is set equal to zero. The second equation consists of 
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the previously defined linear combination of the coefficients represent- 
( ed by cy and T, and is also set equal to zero. The two resulting 
rs simultaneous ordinary differential equations are then expanded and the 
= coefficients of the derivatives regrouped so that equations (2-37) and 


(2-38) can be conveniently expressed in matrix format as follows, 


- p4¢ 03G D26 0G 
ml a maf wf 
( o*H o3H 02H DH 


As it turns out, the equations are a pair of coupled homogeneous differ- 
i ential equations that can be solved as an eigenvalue problem. This 
technique is discussed in the chapter on Numerical Methods. 

The matrices which appear in equation (2-40) are 2 X 2 matrices and 
are summarized in detail in Appendix A. Equation (2-40) and the re- 
spective matrix elements contained in Appendix A are the generalized 
vorticity transport equations for al] angular wave numbers n, at all 
axial wave numbers a, and al] Reynolds numbers. As will be explained in 
detail in the chapter on Boundary Conditions, the generalized equations 


are transformed by appropriate changes of variable depending on the 


angular wave number. The changes of variables from the eigenfunctions 
; G(r) and H(r) to the new eigenfunctions P(r) and Q(r), respectively, are 
determined by the original boundary conditions that must be satisfied, 
( 


particularly at the pipe axis. 
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D0. EQUATIONS FOR n = 0 


For the case n = 0, equation (2-40) is greatly simplified in that 
the two equations uncouple and allow an independent investigation of 
either eigenfunction, G(r) or H(r). Miscellaneous trial calculations 
made earlier suggest that the solutions for G(r) are probably stable for 
all values of a and Re. Therefore, only the solutions for H(r) were 
investigated. It should also be noted that the examination of the axial 
perturbation velocity is of particular interest in this analysis and, in 
the case for n = 0, is not a function of G(r). The change of variable 


required to satisfy the boundary conditions is, 
H(r) = r Q(r) (2-41) 


Taking the derivatives of H(r) yields, 


OH =r DQ +Q 

D2H = r 07Q + 2 0Q 

D?H = r 03Q + 3 02Q 

O*H = r D4Q + 4 03Q (2-42) 


With the substitution of equations (2-41) and (2-42) into equation 


(2-40), the following expression is obtained, 


Mg D*Q + Ms 03Q + Mz 02Q + My DQ + Mo Q = y (Ng 07Q + Nz DQ + No Q) 


(2-43) 
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The primed coefficients are actually the elements from row 2 and column 


2 of the newly transformed coefficient matrices and are not enclosed in 
matrix brackets. However, in general, the new coefficient matrices [M7] 
and [Nj], (where the indices i = 0, 1, 2, 3, 4 and j = 0, 1, 2), are 
formed after the change of variable is made and the terms are collected. 
The primed coefficients used to solve the transformed vorticity trans- 
port equation in terms of Q(r) for n = 0 are defined in Appendix 8, Part 
A. 


—. EQUATIONS for n 


i 
~ 


For the case n = 1, eigenfunctions G(r) and H({r) do not become un- 
coupled and must be solved for in a system of simultaneous equations in 
a form similar to equation (2-40). An additional parameter H(0) is 
introduced into the equations by the change of variables. This parti- 
cular change of variables is required because of the complicated nature 
of the boundary conditions at the axis and is thoroughly explained in 
the chapter on Boundary Conditions. For n = 1, the change of variables 


required to satisfy the boundary conditions are, 


G(r) = -i H(Q) + r2 P(r) 


(2-44) 
H(r) = H(O) + r2 Q(r) 
Taking the derivatives of G(r) yields, 
0G = r2 OP + 2r P 
(2-45) 


02g = r2 2P + 4r DP + 2 P 
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and taking the derivatives of H(r) yields, 


DH = r2 DQ + 2rQ 
02H = r2 02Q + 4r 09 +2Q 
O3H = r2 03Q + Gr D2Q + 6 0Q 
D*H = r2 0*Q + Br D3Q + 12 02Q 


Remember that in general D(H(0)) # DH(O), where OH(0) is the value of 
the first derivative of H(r) evaluated at the point r = 0. Since 
D(H(O)) = 0 here, it vanishes from the expressions for the derivatives 
of G(r) and H(r). In order to accommodate H(0) in the system of equa- 
tions where it appears explicitly, two additional 2 X 1 column matrices 
are required, namely [Ms] and [N,;]. After the changes of variable have 


been made, equation (2-40) now becomes, 


D2p op p 
+ [Ms] {H(0)} = y (Gs \. CNy] | + [No] + (N3] (4091) 
02Q DQ Q 


(2-47) 


The coefficient matrices of equation (2-47) used to solve the trans- 
formed vorticity transport equations in terms of P(r) and Q(r) for n=1 


are defined in Appendix B, Part B. 


F. EQUATIONS FOR n = 6 


For the case n = 6, eigenfunctions G(r) and H(r) remain coupled. 
The change of variables does not introduce any additional parameters and 
the system of simultaneous equations is solved in the form of equation 
(2-40) once again. The change of variables required to satisfy the 
boundary conditions for n = 6 and additionally for all angular wave 


numbers n> 6 are, 


G(r) = r* P(r) 
H(r) = r3 Q(r) (2-48) 


Taking the derivatives of G(r) yields, 


DG = r* DP + 4r3 Pp 
D2G = r* D2p + Sr3 OP + 12r2 P (2-29) 


and taking the derivatives of H(r) yields, 


DH = r? 0Q + 3r2 Q 

02H = r? 02Q + 6r2 DQ + Gr Q 

07H = r3 03Q + 9r2 02Q + 18r DQ + 6 Q 

O*H = r? 0*Q + 12r2 03Q + 36r 02Q + 24 DQ (2-50) 


Substituting equations (2-48), (2-49), and (2-50) into equation (2-40) 
results in the transformed coefficient matrices for the system of 
equations. The coefficient matrices of the functions P(r) and Q(r) and 
their respective derivatives for n = 6 are defined in Appendix 8B, 


Part C. 
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\ III. BOUNDARY CONDITIONS 

S The incorrect formulation of the boundary conditions for the linear- 


re ized vorticity transport equations at the pipe axis has probably been 
the primary reason that previous investigators have failed to predict 
the onset of flow instabilities in circular pipes at any Reynolds number. 
Gawain (Ref. 13] suggested that if these boundary conditions are formu- 
lated in a rigorous and systematic fashion and applied to the vorticity 


eer: Fe a © # @ 
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transport equations, a correct numerical solution may be computed which 
agrees with experimental results. It can be seen from equations (2-40) 
and from the coefficient matrices of Appendix A that the two vorticity 
transport equations are coupled, except for the case n= 0. The re- 
sulting pair of differential equations are second order in G(r) and 
fourth order in H({r). In order to obtain a determinate solution of the 
equations, a total of six boundary conditions are required. More 


specifically, at least two of the boundary conditions must involve G(r) 


> or its derivatives and at least four of the boundary conditions must 


2 


involve H(r) or its derivatives. It might be expected that when the 
equations are coupled, the boundary conditions may also be coupled. 
That is to say that the boundary conditions may consist of various 
linear combinations of the functions and their derivatives, which turns 


out to be the case and is shown later. 
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The boundary conditions at the pipe wall are determined in a fairly 
straight forward manner. They are derived from the fact that the com- 
ponents of the perturbation velocity must vanish at the wall. From the 


definition of the perturbation velocity, defined previously as, 
vrux® (2-27) 


it can be shown that the following three boundary conditions result. 


G(1) =0 
H(1) =0 
DH(1) = 0 (3-1) 


Since six boundary conditions are required and three have been deter- 

mined at the pipe wall, there remain three boundary conditions to be 

determined at the pipe axis. Of these three conditions, at least one 
ie must involve G(r) or its derivatives and at least two must involve H(r) 
or its derivatives. 

The proper derivation of the three boundary conditions at the axis 
turns out to be a non-trival problem. As can be seen in Table 2-1, the 
highest negative power of r that appears in the equations is r* 
Closer inspection of the vorticity transport equations of equation 
(2-40) and the coefficient matrices in Appendix A reveals that the 
equations contain terms in rn, aa e, r, and oe At first glance, 
it may appear that these equations are not satisfied in the limit as ° 
approaches zero. It should be noted that it would be incorrect to 


multiply the equations through by any positive power of r since that 


’ would alter the degree of the singularity of the equations at the axis. 
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This problem can be resolved and subsequently the boundary conditions at 
the axis can be deduced rigorously by expanding the functions G(r) and 
H(r) as power series in r. The highest negative power of r is r* and 
appears in the coefficient matrix (M)] for the pair of functions G(r) 
and H(r). The series expansion for G(r), and similarly for H(r), is 
carried to the fourth derivative as follows, 


G(r) = G(0) + DG(O)r + 2G(0) Sr + D°G(0) Fr + DAG(O) Ey + ..... (32) 


It can be seen that when the series expansion expressions for G(r) and 
H(r) are substituted into equation (2-40), the higher order terms of 
equation (3-2) that contain powers of r greater than four will vanish in 
the limit as r approaches zero. For the first derivatives of G(r) and 
H(r), the highest negative power of r that appears in the coefficient 
matrix (M,) is r. Therefore, the series expansion for the first 
derivative of G(r), and similarly for H{r), up to the fourth derivative 


term is, 


DG(r) = 0 + DG(O) + D2G(0)r + 05G(0) Fr + O4G(0) H+ ..... (3-3) 


The remaining derivatives of the functions are expanded in a similar 
manner through the fourth derivative term. The resulting power series 
approximations for G(r) and H{r) and their derivatives are then substi- 
tuted into equation (2-40). The coefficients are then regrouped in 


ascending powers of r. 
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As mentioned previously, the equations must be satisfied in the 
init as r approaches zero. Therefore, the coefficients of the func- 
tions that contain rae aa cae -. and r must all be zero. The 
remaining coefficients will contain positive powers of r and will vanish 
in the limit as r approaches zero. By considering only the terms con- 
taining negative powers of r and setting those coefficients equal to 
zero, five pair of equations in ten unknowns remain so that the boundary 
conditions can be determined. The unknowns are the functions of G(r) 
and H(r) and their first four derivatives evaluated at the pipe axis, 
r= 0. It is convenient to summarize these equations in matrix format 


below, 


: G(0)) (0 
i oa] =| | on 


H(0) 0 


(3-5) 


a ae ; E 


Re ©?" Lona} lo 


: 026(0)) (0 G(o)) (0 
(r=) [C;] | | = Jol + CialC,] - y [D4]) hae = | | (3-6) 
o2H(0)) lof luo) lo 


a 03G(0)} ia 0G(0)) 0} 
meen) (Cs) +n (aia (1 = =) : :) [Cg] hen a8 (3-7) 
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D*H(0) 02H(0) 


ae 0*6(0) 02G(0) 
(=) tc7] + (ia [Ce] - y [081) 


j (0) 
+ (ia [Cy] - ya? [09]) a -‘ | (3-8) 


The 2 X 2 matrices which appear in equations (3-4) through (3-8) are de- 
fined in detail in Appendix C. Special conditions at the axis arise 
when the determinants of these matrices are evaluated at specific 
angular wave numbers. They are summarized in Appendix 0. 

For each specific value of angular wave number n, it is possible to 
derive a set of boundary conditions at the axis. Special conditions at 
the axis occur for angular wave numbers that cause the coefficients in 
: equations (3-4) through (3-8) to equal zero or the determinants in 
i Appendix D to be zero. That is to say, if a coefficient vanishes or the 
pair of equations are linearly dependent, the functions may be arbi- 
trarily specified and are in general not equal to zero. But if the 
= coefficient, appearing in equations (3-4) through (3-8) is non-zero and 
the determinant exists, the pair of functions must be identically zero. 
This situation occurs for all five pairs of variables when the angular 
wave number n 2 6. 

By taking these special conditions into account for the case n = 0, 


it can be deduced from equations (3-6) and (3-8) that, 


(G¢0)) 2 )02G(0) } _ (9) 


[Hcoyf — [ozH¢0f of or 
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Carrying out this procedure for n = 1, 2, 3, 4, 5, and 6 yields a sim 
plified set of boundary conditions at the axis. The following sets of 
equations (3-10) through (3-16) are essentially the same as those of 
Gawain [Ref. 13] but were independently checked and appear here with a 


few minor corrections. 


n=0 
G(0) = H(0) = 0 a 
a (3-10) 
02G(0) = 0 D2H(0) = 0 (b) 
n=] 
G(0) + i H(0) = 0 
(a) 
0G(0) = DH(0) = 
D3G(0) = 0 D3H(0) = 0 
(3-11) 
e- (ty need (ia [Cg] - y [D9]) da (b) 
- + (j - 
Re losucoyf Deno) 
G(0) 0) 
+ (ia [Co] - y a? [D9]) = 
H(0)} of 
n=2 
G(0) = H(0) = 0 
(a) 
DG(0) + i DH(O) = 0 
wan neon a nnn na nnn no nnn n nn nn nnn nnn nnn e enn nee n ene en nee (3-12) 
02G(0) = 0 D2H(0) = 0 ie 
0*G(0) = 0 D*H(0) = 0 


4 eae Sen ken eee eee ee ee ee ee ewe ei ery Saw gw GMa ene me ese co Fre ee ee ee tee ete eee, 
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a 
a n= 3 
2 G(0) = 0 : H(0) = 
ome (a) 
. 0G(0) = 0 DH(0) = 0 
FF atten mnenanennnnnenennnnnneenananawenneanenenneceneneeenen (3-13) 
ra 02G(0) + 1 02H(0) = 0 ‘ee 
S 03G(0) = 0 D3H(0) = 0 
- n=4 
. G(0) = 0 HCO) = 0 
: 0G(0) = 0 DH(O) = 0 (a) 
02G(0) = 0 02H(0) = 0 
ee eee (3-14) 
a 0?G(0) + i D3H(0) = 0 ee 
0*G(0) = 0 0*H(0) = 0 
: n= 5 
" G(0) = 0 H(0) = 0 
0G(0) = 0 DH(O) = 0 (a) 
026(0) = 0 D2H(0) = 0 
serene enone 2-22-22 22-2 2-22 o-oo no-no (3-15) 
D3G(0) = 0 -  D3H(0) = 0 ™ 
0*G(0) + i D*H(0) = 0 
n26 
5 G(0) = 0 
- DG(0) = 0 H(0) = 0 
: D2G(0) = 0 DH(O) = 0 (a) 
7: 03G(0) = 0 02H(0) = 0 
ween enn nnn n nnn ne nn nnn enn nnn nnn nnn onan nnn one n nnn n-ne eee (3-16) 
0*G(0) = 0 D3H(0) = 0 
ne (b) 
4 0*H(0) = 0 
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Recall that exactly three boundary conditions are required at the 


axis to finalize the solution of the vorticity transport equations. 
Equations (3-10) through (3-16) contain from four to ten boundary con- 
straints including linear combinations of several boundary conditions, 
as suggested earlier. By introducing appropriate changes of variables, 
the pair of functions in G(r) and H(r) and their respective derivatives 
can be transformed to a new pair of functions in terms of P(r) and Q(r) 
and their derivatives. The transformed system of equations take the 
form of equation (2-40) for all values of n, except for n= 1 where 
equation (2-47) is applicable. The change of variables is dependent on 
each specific value of n up through six and is chosen in such a way that 
the boundary conditions in subset (a) of equation (3-10) through (3-16) 
are satisfied identically. The boundary conditions in subset (b) must 
be satisfied explicitly in the finite difference equations near the 
axis. Development of the finite difference equations is discussed in 
the chapter on Numerical Methods. 

The following sets of equations are taken from Gawain [Ref. 13] and 


are included here so that a complete treatment of the subject on bound- 


ary conditions is contained in this paper. 


_ er 2 5 R ’ =: Fa uy : a ced — | 5 News Tae oe 2 
Tae a aT at a he a ee eae ee ete ele Ee eee ee tee 
PSs 


. 
n=Q 
Change of Variables 
G(r) = rP(r) 
(a) 
ae H(r) = rQ(r) 
F- Boundary Conditions at Axis (3-17) 
: DP(O) = 0 09(0) = 0 
: (b) 
07Q(0) = 0 
x Boundary Conditions at Wall 
: P(1) = 0 Q(1) = 0 
(c) 
0Q(1) = 0 
The solutions for P(r) and Q(r) become uncoupled for n = 0 
n=l 
! Change of Variables 
, G(r) = -1 H(O) + r2 P(r) es 
‘ a 
H(r) = H(O) + r2 Q(r) 
Boundary Conditions at Axis (3-18) 
OP(O) = 0 0Q(0) = 0 
56 02°(0) (Pco)) {-i 
> x [C7] + 2ia(Cg] + ja[{Cg] H(0) 
| (n2Q¢0) (q¢0) (a 


Boundary Conditions at Wall 


- 7 H(0) + P(1) = 0 H(O) + Q(1) = 0 
(c) 
-2H(0) + 0Q(1) = 0 


= 


Ne 


r, 
Jae 8 


~s  e  D ahaa”) 7 
iv. eee, eee Seer Tana ee 
. 


x 
s nz2 
i Change of Variables ; 
= G(r) = - { r DH(O) + r2 P(r) 
(a) 
H(r) = r DH(O) + r2 cr) 
: Boundary Conditions at axis (3-19) 
. P(0) = 0 Q(0) = 0 
i (b) 
. D2P(0) = 0 02Q(0) = 0 
: Boundary Conditions at Wal] 
- 4 DH(O) + P(1) =0 DH(O) + Q(1) = 
a: (c)- 
-DH(O) + 0Q(1) = 0 
n= 3 
: Change of Variables 
i G(r) = r2 P(r) 
G (a) 
: H(r) = r2 Q(r) 
Boundary Conditions at Axis (3-20) 
P(O) + 1 Q(O) = 0 
(b) 
OP(O) = 0 0Q(0) = 0 
Boundary Conditions at Wal} 
P(1) = 0 Q(1) = 
4 (c) 
0Q(1) = 
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( n2z4 : 
. Change of Variables 
: G(r) = r3 P(r) 
xg (a) 
: H(r) = r Q(r) 
re Boundary Conditfons at Axis (3-21) 
: P(O) + 1 Q(0) = 0 
s (b) 
7% OP(O) = 0 0Q(0) = 0 
“ Boundary Conditions at Wall 
e P(1) = 0 Q(1) = 0 
ba (c) 
i 0Q(1) = 
s n=5 
hg Change of Variables 
ss G(r) = r3 P(r) 
. (a) 
H(r) = r3 Q(r) 
Boundary Conditions at Axis (3-22) 
P(O) = 0 Q(o) = 0 
(0) 
3 DP(0) + i 0Q(0) = 0 
- Boundary Conditions at Wall 
Hi P(1) = 0 Q(1) = 0 
es (c) 
< 0Q(1) = 0 
¢ 
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n 26 
Change of Variables 


G(r) = r* P(r) 


(a) 
H(r) = r? Q(r) 
Boundary Conditions at Axis (3-23) 
P(0) = 0 Q(0) = 0 
(b) 
PALATE CRE ee 0Q(0) = 0 
Boundary Conditions at Wal} 
P(1) = 0 Q(1) = 0 
(c) 
0Q(1) = 0 


It should be noted that the change of variable relations in subset (a) 
of equations (3-17) through (3-23) does not change the order of the 
derivatives appearing in the transformed vorticity transport equations. 
Upon close examination of subsets (b) and (c) of these equations, it is 
apparent that the correct number of boundary conditions required for the 
solution of the transformed equations has resulted from the change of 
variables. The vorticity transport equations are now expressed in terms 
of P(r) and Q(r) and their derivatives and have exactly three boundary 
conditions at the pipe wall for all values of n and three boundary 
conditions at the axis for all values of n, except for n = 1 and 2. The 
fact that there are four boundary conditions at the axis for angular 


wave numbers n = 1 and 2 warrants further comment here. 
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Notice in equations (3-18a) and (3-19a) for the cases n = 1 and 2 
that an additional parameter is introduced by the change of variable 
equations. The additional parameter for n = 1 is H(0) and for n = 2 is 

; DH(O). The introduction of the parameter for each case requires the 
specification of four boundary conditions at the axis instead of the 

. usual three. This conclusion was reached by Gawain [Ref. 13] and is 

analyzed in detail in his paper. There is one other peculiarity about 
the boundary conditions at the axis for n = 1 that should be mentioned. 
The axis boundary conditions in equation (3-18b) contain a pair of 
equations that involve the complex perturbation frequency y, which 
represents an eigenvalue of the solution for the equations. The pres- 
ence of y in the axis boundary condition for n = 1 is an important 
factor in formulating a solution for the vorticity transport equations 
and is included in the eigensystem of the finite difference equations. 
The details of the implementation are discussed later. 

The solution of the problem for predicting the onset of flow in- 
stabilities in circular pipes is now at hand. The original vorticity 
transport equations (2-40) are transformed to the functions P(r) and 
Q(r) and their respective derivatives by the appropriate change of 
variables depending on the specific angular wave number being investi- 
gated. After imposing the rigorously deduced boundary conditions, a 


solution can be determined in the perturbation quantities. 


7 5 a : — 
mE a Se 


IV. PERTURBATION VELOCITY 


The solution of the transformed vorticity transport equations de- 
veloped in the previous chapter can now be expressed in terms of the 
perturbation quantities P(r) and Q(r). The original functions G(r) and 
H(r) are actually components of the vector potential function W defined 
in equation (2-28), 


W = (F(r) é, + G(r) é. + H(r) 1 e* (2-28) 


The axial component F(r) was previously set to zero since one function 
could be arbitrarily specified. From equation (2-27) the perturbation 


velocity vector was expressed as, 
=v xW (2-27) 


Therefore the functions G(r) and H(r) represent the behavior of the 
perturbation velocity in the flow field. Referring back to the defini- 
tion of the curl operation in equation (2-39), the relation for V x W 
has already been set up. The resulting expression for the perturbation 


velocity is given by, 


vated Sasa te Me dee tote te eee TT Te an ee te he Seng Ae Bete ne Sas Ye ye peer ns Te 
. 
. 


i 
; The components of the perturbation velocity are, 

§ u = DH(r) + ry - in G(r) (4-2) 
: v = fa H(r) (4-3) 
§ w= fa G(r) (4-4) 
: Digressing for a moment, it can be readily seen that since the pertur- 
3 bation velocity necessarily vanishes at the pipe wall, r = 1, the 
- boundary conditions at the wall in equation (3-1) can be deduced. 

3 Recall that the functions P(r) and Q(r) appear as the result of the 
: change of variables for each specific angular wave number contained in 
subset (a) of equations (3-17) through (3-23). Analysis of the axial 
. component of the perturbation velocity u becomes useful in interpreting 
s the behavior of the functions P(r) and Q(r) along the radius of the 


pipe. Equations for the axial perturbation velocity in terms of P(r) 
and Q(r) are derived from equations (4-2) through (4-4) with the change 
of variable applied for each angular wave number being investigated and 


appear as follows, 


: n=0 u=2Q+r 0Q (4-5) 
; n=l: u=3rQ+ter20Q-irP (4-6) 
n=6: u=4r2Q + r30Q - 16 r3 P (4-7) 


The axial perturbation velocity is computed in part VII of the main 
investigative program for each angular wave number. To prepare the 


perturbation velocity u for plotting versus pipe radius, it is first 
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normalized. Since this velocity is in general a complex quantity, the 
re perturbation velocity with the largest’ magnitude is determined by, 

IU. = (up? + uy2)" | (4-8) 
MAX 

Pe The normalizing velocity constant will then become U. and the maximum 
2 normalized velocity will become, 
( 
= uz 1+ 1(0) (4-9) 
H u 
a c 
The remaining normalized velocities are found by dividing through by Un. 
Ce 
a The normalized perturbation velocity is then plotted versus the normal- 
: ized pipe radius for various axial wave numbers and Reynolds numbers at 
i the three angular wave numbers investigated. The plots and results are 
4 
i discussed in detail in the chapter on Results. 
4 . 
é 
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V. NUMERICAL METHODS 


A. GENERAL METHODS USED 


In general, the vorticity transport equations are a pair of coupled 
homogeneous differential equations which are second order in the pertur- 
bation quantity G(r) and fourth order in the perturbation quantity H(r). 
A change of variables is introduced so that the number of boundary 
conditions at the axis is reduced to three for each angular wave number 
investigated, except for n = 1. Recall for this case that there are 
four boundary conditions imposed at the axis because of the introduction 
of the additional parameter H(0). By applying the change of variables 
to the original pair of equations, the vorticity transport equations are 
transformed into an equivalent pair of coupled homogeneous differential 
equations that are second order in P(r) and fourth order in Q(r). The 
solution of the vorticity transport equations will now be in terms of 
the perturbation quantities P(r) and Q(r) which are also referred to as 
eigenfunctions, and y which is an eigenvalue. The axial perturbation 


velocity is expressed in terms of P(r) and Q(r) as derived previously 


and is examined to determine the behavior of the perturbation quantities. 


The resulting eigenvalue y is of primary interest since the magnitude 
and sign of the real part of y determines the relative stability or 
instability of the pipe flow problem. 

A numerical solution of the transformed vorticity transport equa- 
tions is now sought. The convenient matrix format of the equations is 


given below. 
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o*p p3P p2p oP } p 
[Mg] + [Ms] + [Mo] + [My] + (Mo) = 
0*q psq 02q 09 Q 
p2p op p 
y (Ng) + (Ny] + [No] (5-1) 
02Q DQ Q 


The details of the transformed couffictent matrices for n=0, 1, and 6 


are given in Appendix B. It is ‘convenient to approximate the deriva- 
tives of P(r) and Q(r) by a finite number of discrete unknowns along a 
radial mesh consisting of N interior points. It is not possible to 
compute the continuous values of the functions and their derivatives 
along the pipe radius using a digital computer. Therefore, discrete 
finite difference techniques are used. For the problem at hand, the 
normalized nondimensional radius of the pipe, which becomes the in- 
dependent variable, is divided into a computational mesh. The standard 
method of representing this mesh is with uniformly spaced stations over 


the interval 0 < r< 1, Figure 5-1. 


Pipe Pipe 
wae Wall 
r=] 


N-3 N-2 N-1l ON 


Figure 5-1. Radial Mesh, Standard Method 
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oy The vorticity transport equation coefficients are then evaluated at each 


r station of the radial mesh resulting in N uniformily spaced discrete 

be values with N + 1 equal intervals and N + 2 total points, including 

: r=0 and r=1. The spacing between the stations is denoted by 
h = 1/(N + 1). 


A refinement of the method for generating a computational mesh is 

referred to as the half station method and is used in this numerical 

e analysis. All the radial stations of Figure 5-1 are moved to mid- 
interval or one half station (h/2) toward the axis with an additional 


station appearing at one half station from the pipe wall, as shown in 


Figure 5-2. 

Pipe Pipe 
.. Axis h/2 h h h/2 wall 
i r=0 i aig r=1 
. 1 2 3 4 #5 N-4 N-3 N-2 N-1l ON 

Figure 5-2. Radial Mesh, Half Station Method 

Therefore for N interior stations, the method results in N intervals and 

N + 2 total points including the boundaries. The first station, r,, is 
> a distance h/2 from the axis and the Nth station, Tye is a distance h/2 
i from the wall. The general] expression for the ith station, correspond- 
he ing to the ith increment of the radius, is given by, 
d 


rozh(i- 3) - 4 =1, 2, 3....N (5-2) 


where h = 1/N. 
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An advantage of the half station method over the standard method is 
that there is better resolution at the stations near the boundaries 


where instabilities are thought to originate and where boundary layer 
effects are a factor. Therefore, the behavior of the perturbation 
quantities can be examined more closely. Additionally, better approxi- 
mations for the derivatives of the functions P(r) and Q(r) can be 
realized with a smaller incremental distance to the points immediately 
adjacent to the boundaries. Even better resolution at the boundaries 
can be achieved by using a nonuniform computational mesh. Arnold 
: {Ref. 11] introduced a change of independent variable by means of hyper- 
z bolic functions which included a mesh offset parameter so that he could 
contro] the concentration of mesh points at either the axis or wall. 
The method would certainly be useful in reducing computational time by 
reducing the total number of mesh points required for a satisfactory 
solution of the governing equations. However, the full implications of 
the nonuniform mesh with respect to a and y dependence on the mesh 
offset parameter have not been fully explored. While this method is not 
used here, it may well prove to be useful in follow-on numerical in- 
vestigations of the pipe flow problem. 

The solution for the vorticity transport equations in terms of the 
perturbation quantities P(r) and Q(r) are now represented by a finite 
number of discrete unknowns. The governing differential equations are 
rewritten in finite difference form for N stations along the mesh re- 
sulting in a finite number of unknowns namely, P,, Po, Pa...Py and Q;, 


Q; Qs...-Qy- In general, this produces 2N simultaneous equations in 2N 


eet ue ef Toay Mowe en reat ny 


unknowns. The details of the derivation of the finite difference equa- 


tions are explained shortly. Substituting the finite difference equa- 


tions into the vorticity transport equation matrix format of equation 


(5-1) will result in N pairs of coupled simultaneous equations. These 


equations are regrouped and subsequently solved in the following eigen- 


value problem format, 


CA] {X} = y (B] {X} (5-3) 


The column vector {X} represents the rearranged unknowns as 


Py 


The symbol y denotes the characteristic value or eigenvalue of the 


solution. Details of the composition of the [A] and (B8] matrices are 


discussed for each specific angular wave number investigated. 
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B. FINITE DIFFERENCE EQUATIONS 


The use of finite differencing techniques is the basis for solving 
the governing vorticity transport equations in the eigensystem equation 
(5-3). It is required that the governing equations must be satisfied at 
each of the interior stations i = 1, 2, 3....N, therefore producing 2N 
equations in 2N unknowns. 

Except for stations adjacent to the boundaries, the derivatives of 
the perturbation quantities in the governing equations at the interior 
stations are usually approximated by standard central finite difference 
equations. These equations were taken from Table 7.27 of Ketter and 


Prawel (Ref. 14] and are summarized below, 


eae n = 1 5 of 


° 


1 f_ : ae Lee ae eee 


2 
O°F. 


1 7 
3 = - - - 477 
O°Fs = ane Fas BF ing * 135.1 ~ MSFi4y * BFing ~ Fini * Tam NO'F; 


(5-7) 
Jig Ne : : = 
OF = oe | Fig * 1A; 2 7 39F;1 + 56; ~ 39F,) + FD - Fiag 
+ she hADsF (5-8) 
zag OF; 


where F is an arbitrary function. 
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As can be seen from the equations above, the derivative of a func- 


tion at a particular station can be approximated by a linear combination 
of the discrete unknown functions about a central point. The term 
outside the brackets in each equation represents the truncation error of 
the corresponding finite difference approximation. Notice that the 
truncation error of equations (5-5) through (5-8) are all of order h* 
(Oh*). It was hypothesized and shown by Gawain and Ball (Ref. 15] that 
finite difference approximations of consistent order truncation error 
reduced the overall error of the numerical solution in a typical problem 
by a factor of five for truncation error Oh?. For the present numerical 
analysis, a computational mesh of N = 50 was used for all cases. Now 
with the interval h = 1/N = 0.02, examination of the consistent fourth 
order truncation error term will yield an error of the order of 5x10"? 
for all four derivatives approximated, since the magnitude of the 
normalized functions and derivatives are of order one. for an interval 
size of h = 0.02, a consistent truncation error Oh* should significantly 
reduce the overall error of the solution. 

Special problems arise in the derivation of the finite difference 
equations at stations near the boundaries. Several approaches can be 
taken to approximate the derivatives at these stations, which include 
backward, forward and central finite differencing techniques. Unfor- 
tunately, the central method previously discussed will involve imaginary 
stations outside the boundaries for stations close to the axis and wall, 
i.e., stations i = -1, -2, -3, and N+l1, N+2, N+3 are involved in the 


fourth derivative approximation at station i = 1 and j = N, respectively. 


Unlike the boundary conditions for many problems well suited for this 
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type of numerical analysis, the vorticity transport equation boundary 
(| conditions for the pipe flow problem are very complicated by comparison. 
The boundary conditions for each angular wave number are different; some 
. involve additional parameters while others contain the eigenvalue y. In 
other words, the conditions at the imaginary stations for conventional 
problems can normally be expressed in terms of conditions at the real 
stations in the computational mesh. It is therefore impractical to use 
the central method for the pipe flow problem because of the complicated 
nature of the boundary conditions. 

The method of forward differencing fis used near the pipe axis and 
backward differencing is used near the wall. These two differencing 
methods will be referred to, collectively, as noncentral finite differ- 
encing techniques. It will be shown later that the two methods are 
essentially equivalent. A rigorous derivation of the noncentral finite 
difference equations is based on the Taylor series (power series) expan- 
sion of the function at a point. Since the governing equations contain 
the fourth order derivative of Q(r), the derivation for this case is 
examined in detail here. 

The Taylor series approximation of Q;(r) at station i = 1 appears as 
equation (5-9). 


h,2 h,3 h 
h (5) (5) (5) 
Qi = Q(0) + 0Q(0) 3+ D2Q(0) ar it D°Q(0) -- + 0*Q(0) 


Z h,5 h,®& h,? 
A (3) (3) (3) 
+ 05Q(0) “ar + 0°Q(0) a | al 07Q(0) > (5-9) 
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The approximation for Q, at stations {i = 1, 2, 3, 4, 5, 6 are expressed 


J . _ in convenient matrix format below, 
(1. GD GB GB ® GB 
rr) ee ie i ie 
; @» & ¢ & 
p< I a 
Q 
5,2 
Q2 1 5 (3) ° 
Qs 2; 0 
Qs 
pa Q 7 Coy e ° 
: t 30 ° 20h 
% Qe 
> 
> 2. ae 
‘ 11,2 
eae. 
2 a 
; Q(0) 
: hDQ(0) 
7 h202Q(0) 
: h303Q(0) 
h40*Q(0) * {£1} n®pegco) (5-10) 
h5D5Q(0) 
h8D8Q(0) 
h7D7Q(0) 


where the dots represent matrix elements that can be easily deduced. 
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Since the fourth order derivative of Q(r) is present in the vorti- 
city transport equations, there will be two boundary conditions speci- 
fied at the axis and two at the wall. It can be shown that two columns 
of the matrix in equation (5-10) can be eliminated since two of the 
functions in the column vector are either zero or can be expressed in 
terms of additional parameters. The matrix, which will be denoted as 
[AA], now becomes a square 6 X 6 matrix. Notice that the column vector 
{E,} represents the last nonvanishing term in the Taylor series. The 
approximation for the first derivative of Q;(r) at station i = 1 is 


given by, 
' Co cy" cb)* 
0Q, = DQ(0) + 07Q(0) 5 + 03Q(0) = | aga 0*Q(0) a" * 0®Q(0) <7 
)° oy 
+ D6Q(0) “eT + D07Q(0) —I- fo hae (5-11) 


The Taylor series expansion for the second, third, and fourth deriv- 
ative can be found by successive deferentiation of equation (5-11). The 
first four derivatives of Q; at station i = 1 can be expressed in matrix 


format as, 
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z , ® @ ¢* @ 6° 
1@ 3 a a : 
pe 1.2 ,1,3 ,1.¢ 1,8 
‘ j@) @& td G& 
7 ie 001 @ Sr ar ar or 
h?02Q, 
eh ib . & @ ¢* 
: h4p*Q, aS ar es ee | oe 
§ 
1,2 ,1,8 
(i) G&) 

: 000 0 1 @) i 4 
: 4x8 
7 Q(0) 
: h0Q(0) 

h202Q(0) 

h?03Q(0) 

h*D4Q(0) + {E2} h8D8Q(0) (5-12) 

h505Q(0) 

h®0®Q(0) 

h707Q(0) 
. It can be shown here as well that two columns of the matrix in 
i equation (5-12) can be eliminated so that the matrix denoted [CC] jis now 


a 4 X 6 matrix. The column vector {E,} also represents the last non- 
vanishing term in the Taylor series. It can be seen that the approxi- 
mations for the derivatives of Q; at stations ji = 2 and i = 3 can be 
determined by replacing (1/2) in the [CC] matrix with (3/2) and (5/2), 


respectively. 
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To show how the noncentral finite difference approximations near the 


axis are formulated, the case for angular wave number n 2 6 is examined. 


The boundary conditions at the axis for the function Q are, 


Q(0) = 0 and 0Q(0) = 0 (5-13) 


Applying these conditions to equations (5-10) and (5-12) will eliminate 
the first two columns of the [AA] and [CC] matrices since the first two 
terms of the column vectors vanish. The truncation error terms {E,} and 
{E.} are very small quantities and are neglected at this point. In 
shorthand notation, equation (5-10) now becomes, 


h202Q(0) 
h303Q(0) 
h40*Q(0) 
h§05Q(0) 
h®D&Q(0) 
h707Q(0) 


(5-14) 


and equation (5-12) becomes, 


h?02Q(0) 
h3D3Q(0) 
h40*Q(0) 
h5D5Q(0) 
h®D&Q(0) 
h?07Q(0) 
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. 
% Equation (5-14) is aultiplied through [AA] = and becones, 
re h202Q(0) Q 
a h?D2Q(0) 1a 
cr: h*0*Q(0) ay Qs 
. = [AA] (5-16) 
h505Q(0), 6x6 Q& 
Pe nensq(o) | Qs 
Z h?07Q(0) Qs 
and is then substituted into equation (5-15) giving, 
‘ 
‘i Qs 
: nO, Qe 
a h2020, a Qs 
iB = [CC] [AA] (5-17) 
h?03Q, 4x6 6x6 G 
h*0*Q, Is 
Qe 
: Now, the fourth derivative of Q; at station i = 1 is approximated in 
: terms of the discrete functions of Q at the interior mesh points with 
truncation error Oh*. The values on the bottom row of the resulting 
: 4X6 matrix after multiplication in equation (5-17) are the coeffi- 
: cients of the column vector Q, through Q, for 0*Q,. In order to have 


consistent truncation error Oh*, the last row and last column of the 
{AA] and [CC] matrices are deleted for the third derivative of Q. 
Matrix [AA] becomes a 5 X 5 and is now inverted and matrix [CC] becomes 


a3X 5. The values on the bottom row of the resulting 3 X 5 matrix are 


the coefficients of Q, through Q, which approximates 02Q,. The rows and 
columns of matrices [AA] and [CC] are successively reduced to determine 
the coefficients for 02Q, in terms of Q, through Q,, and 0Q, in terus of 
Q; through Q, so that consistent truncation error Oh* is preserved 
throughout. The noncentral finite difference equation for 0*Q, will now 
appear as, 


0*Q: = fr (ArQs + AaQe + Aas + AqQs + AsQs + Ase) + Oh = (5-18) 


where A, through Ag are the coefficients taken from the bottom row of 
the resulting matrix in equation (5-17) 

In a similar manner, by imposing the required boundary conditions, 
this process can be carried out to determine the noncentral finite 
difference equations for the first two derivatives of P(r) at several 
points near the boundaries as required. 

It is important to note here that the noncentral finite difference 
approximations near the wall can be formulated by using the same tech- 
niques as those used near the axis. Once again taking the case for 
angular wave number n 2 6, the boundary conditions at the pipe wal! for 
the function Q are, 


Q(1) = 0 and 0Q(1) = 0 (5-19) 


Rather than set up complete new [AA] and [CC] matrices at different 
values of r in the Taylor series expansion near the wall, assume that 


the boundary conditions in equation (5-19) are imposed at the axis, 
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r=0. Notice that the boundary conditions in equation (5-19) are now 


identical to the boundary conditions in equation (5-13). The method for 
computing the coefficients of the noncentral finite difference approxi- 
mations at the wall can be carried out as discussed previously and, in 
this case, have already been computed. It can be shown that the A; 
coefficients of equation (5-18) are applicable to 0*Q, at the wall but 
will appear in reverse order. The noncentral finite difference equation 


for 0*Q, is given by, 


1 
49, = 
O°Qy = AF CAGQy-5 + A5Qy-g * AgQy-g * AgQy-o + AnQy-g + AyQy) + OWS 


(5-20) 


The method for determining the coefficients for the remaining deriv- 
atives of Q and P near the wall proceed as though the boundary condi- 
tions were imposed at the axis; the coefficients are computed, their 
order reversed and coefficients of the odd order derivatives undergo a 
sign change. Order reversal and sign changes are accomplished in part 
II of the main investigative program when the precomputed noncentral 
finite difference coefficients are read in as data. 

As can be seen from the central finite difference approximations of 
equations (5-5) through (5-8), the third and fourth derivatives of a 
function require three points either side of the central point being 
approximated. The first and second derivatives require only two points 
either side of the central point. Therefore the central finite dif- 


ference approximations are applicable only to interior stations i = 4 


through N-3 for 0*Q, and 0°Q;, and stations i = 3 through N-2 for 07Q,, 


0Q;, 02P,, and OP,. This leaves either two or three points immediately 
adjacent to the boundaries to be approximated by the noncentral finite 
difference equations just developed. 

Since the central finite difference approximations are applicable to 
the above mentioned interior mesh points and are valid for all the de- 
rivatives of P(r) and Q(r) as well, the coefficients are computed only 
once in part I of the main investigative program in double precision 
format. The coefficients of the noncentral finite difference approxi- 
mations, which are dependent on the boundary conditions for each angular 
wave number, are computed in an independent computer program. The 
resulting noncentral coefficients for the derivatives of P(r) and Q(r) 
at mesh points near the boundaries are also computed in double precision 
format and are read into part II of the main investigative program as 


data. 
C. SPECIFIC METHODS FOR n = 0, 1, AND 6 


The vorticity transport equations have now been expanded into a set 
of discrete simultaneous equations with respect to a half station 
computational mesh of N points (stations). The derivatives of the 


governing equations were approximated by central and noncentral finite 


difference coefficients, with consistent fourth order truncation error, 


in terms of the discrete unknown functions of the system. The con- 
venient matrix format of the eigensystem in equation (5-3) can now be 
solved on a digital computer using a specialized Fortran subroutine, 
which will be discussed later. It is important to understand the com- 


position of the [A] and [8] matrices and the column vector of unknowns 


{X}. The problem setup is examined in detail for each angular wave 
number investigated. In the following discussions, the vorticity trans- 
port equations are referred to in the transformed state after the change 
of variables to P(r) and Q(r) have been made. The boundary conditions 
in equations (3-17), (3-18), and (3-23), and the coefficients in 
Appendix B also apply. 

For the case n = 0, the pipe flow problem becomes greatly simplified 
in that the equations become uncoupled and only the solution for Q(r) is 


investigated for reasons explained earlier. The governing equation in 


terms of Q(r) is represented by equation (2-38). At the outset, it was 


determined that a mesh size of N = 50 was sufficient to produce a satis- 
factory numerical solution of the problem. It was necessary to deter- 
mine when the solution was relatively independent of mesh size. Prelim- 
inary trial computations were performed for various mesh sizes, N = 25 
to 100. From N = 45 to 100, the solution of the equation for n = 0 did 
not vary appreciably. Additionally, the solution was relatively in- 
dependent of Reynolds numbers below 15,000 for mesh sizes of about 50. 
In order to reduce computation time and the eigensystem matrix size, a 
mesh of N = 50 was used in this analysis. 

In general for angular wave number n = 0, the [A] and (B] matrices 
are 50 X 50 and appear as banded diagonal matrices. For a consistent 
order truncation error Oh* used here, the derivative of highest order 
requires the greatest band width. From the governing equation in Q, it 
can be seen that the highest order derivative in [A] is 0*Q and in (B] 
is D?Q. At stations near the boundaries, the noncentral finite differ- 


ence coefficients require a band width of six in matrix [A] and four in 
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matrix [B]. At interior stations, the central finite difference co- 
efficients require a band width of seven for [A] and five for [(B]. 
Therefore the [A] and [B] matrices of equation (5-3) contain nonzero 
elements as shown in Figures 5-3 and 5-4. The X elements contain non- 
central finite difference coefficients and the 0 elements contain 
central difference coefficients. The X and @ symbols denote the station 
at which the approximation is made. The composition of the column 


vector {X} can also be seen in Figures 5-3 and 5-4. 
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Figure 5-3. [A] Matrix for n = 0 
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Figure 5-4. [B] Matrix for n = 0 


The actual values of the elements of the [A] and [B] matrices are 
computed in part IV of the main investigative program. The calculations 
are rather lengthy and tedious and are accomplished in iterative loops 
in part IV. To save computation time, calculations are made only for 
the ee elements of the matrices. Recall that the governing equa- 
tion in Q(r) for n = 0 was given in equation (2-43). In discrete form 


for the ith station, equation (2-43) appears as, 
Ma; D*Q, + Ma; D°Q; + M2, 07Q; + Mi; 0Q; + Mo; Q; = 
Y(Ng; D2Q, + Ni; 0Q; + No; Q;) (5-21) 
where the primed quantities are element (2,2) of the corresponding 
transformed vorticity transport equation matrices given in Appendix 8, 


Part A. At station i = 1, the top row of the [A] matrix will have six 


elements and the top row of [B] will have four elements. 
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For example, the matrix elements for the first row of [A] can be 

found as follows, where the a, coefficients are analogous to the co- 
4 

efficients A; of equation (5-18) and are equal to a; = A,/h . Statlarly 

b, = B,/h3, Cc, = C,/h?, and d, = D,/h, where by, Cy» and d; are the 

noncentral finite difference coefficients of the remaining three 


derivatives. 


Ara {Qi} = (Ma a, + Ms. b; + Ma. C, + Ma d, + Mo) {Qi} (5-22) 


Ar 2 {02} = (Ma a2 + Ma ba + Mz co + Mi de) {Qe} (5-23) 
Ars {Qs} = (MZ ag + Ms ba + M3 ca + Mi ds) {Qs} (5-24) | 
Ar 4 {Qe} = (MG ag + M5 be + M3 ca) {04} (5-25) 
Ars {Qs} = (Ma as + Ms bs) {Qs} (5-26) 
Ar {Qe} = (Mi a6) {Qe} (5-27) 


where A; , is the matrix element in the first row and the first column 
oa of [A], and so forth. This procedure is carried out for N equations 
along the computational mesh until] all the nonzero matrix elements of 


[A] and [B] are computed. 


E Derivation of the noncentral finite difference coefficients for 
= other angular wave numbers is fairly straightforward when accomplished 
Y by the methods outlined in the previous section. The boundary condi- 
‘ tions for Q at the axis for n = 0 are, 
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. 0Q(0) = 0 and 07Q(0) = 0 (5-28) 


After imposing these conditions, the second and fourth columns of the 
matrices in equations (5-10) and (5-12) can be eliminated since the 
corresponding terms in the Selman: vecters vanish. The coefficients are 
then computed as described before in an independent computer program and 
then read into the main investigative program as data. Notice that the 
bondary conditions for the function Q at the wall are the same for n = 0 
and n 2 6. Therefore, the noncentral finite difference coefficients 
need to be computed only once. 

The case for angular wave number n 2 6 will be examined next since 
it more nearly represents the general case. Since the vorticity trans- 
port equations remain coupled for n 2 6 and no additional parameters are 
introduced, a system of 2N equations in 2N unknowns will result. For a 
mesh size of 50, the [A] and [8] matrices become 100 X 100 and are 
solved in the form of equation (5-3). The first N equations of the 
system are represented by the discrete form of equation (2-37), which 
after transformation by the appropriate change of variables corresponds 
to the top equation of (5-1). The second N set of equations are repre- 
sented by the discrete form of equation (2-38), which after transfor- 
mation by the appropriate change of variables corresponds to the bottom 
equation of (5-1). As a matter of convenience for programming, the [A] 
and [8] matrices are partitioned into four 50 X 50 quadrants, which are 
individually banded diagonal matrices. The properties of the [A] and 


(B] matrices for n 2 6 are summarized in Table 5-1. 
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Table 5-1. Properties of the [A] and [8] Matrices for n 2 6 


HIGHEST ORDER | NONCENTRAL CENTRAL 
DERIVATIVE COEFFICIENTS COEFFICIENTS 
PRESENT REQUIRED REQUIRED 


* Number of stations near the boundaries where the noncentral finite 
difference approximations apply. 


It can be deduced from the data given in Table 5-1 that the (A] and 
[B] matrices will contain nonzero elements as shown in Figures 5-5 and 
5-6. The column vector {X} of unknowns is also given. 

The method for determining the noncentral finite difference co- 
efficients for n 2 6 was discussed in a previous example. 

The case for n = 1 is unique in that an additional parameter H(0) is 
introduced by the change of variables and a pair of coupled boundary 
conditions at the axis in equation (3-18b) contain the eigenvalue y. 
These peculiarities are significant factors in correctly setting up the 
problem. Recall that in the coupled governing equations (2-47) the 


additional parameter H(0) appears explicitly and two additional column 
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XXXXX 1 XXXXXX Q: 
XXXXX III + XXXKXK IV Qs 
00000 1 XXXXXX Qs 
, 00000 + 0000000 
00000 : 
00000 0000000 Q 
0000 | 0000000 ma 
o | 0000000 : 
00000 | XXXXXX Q 
XXXXX J KXXXXX Nel 
XXXXX | XXXXXX QW 
Figure 5-5. [A] Matrix for n 2 6 
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Figure 5-6. [B] Matrix for n 2 6 
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‘ matrices are required in the system of equations. This was due to the 


4 change of variables introduced by equations (2-44). The appearance of 
the eigenvalue y in the pair of boundary equations at the ante comes 
: about from the rigorous method in which the boundary conditions were 
bs derived in Chapter III. The coupled axis boundary conditions in equa- 
é tion (3-18b) are repeated here, 

: 

96 02P(0) P(0) “i 

~ - Fe [Cr] + 2ia(Cg) + fa[Co] H(0) 

a 02Q(0) Q(0) 1 

: (5-29) 


4 . P(0) -j 
. *y (: [Dg] + a2 [D9] Heo} 
* Q(0) 1 


where the coefficient matrices are contained in Appendix C. Notice that 


the parameter H(0O) appears explicitly in these two equations as well. 


Sees Le, 


Additionally, the parameters D2P(0), D2Q(0), P(0), and Q(0) are also 
introduced. Obviously, this poses some problems since there are only 
two additional equations and five unknown parameters. H(0) can be used 


since it appears in the 2N discrete equations of the system. Q(0) is 


in 


chosen as the second variable as a matter of convenience. The remaining 
three functions can be approximated in terms of discrete functions of P; 
in the case for D?P(0) and P(0), and Q; and Q(0) in the case of 02Q(0). 
By referring to the matrix form of the Taylor series expansion of a 
function at a point in equation (5-10), it can be seen how the approxi- 


mations of these functions are made. For the present situation where 
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DP(0) = 0, eliminate the second column of the aatrix and equation (5-10) 
will become, 


Py P(0) 

1 Pe h202P(0) 

Ps? = CAA} _-< h3D3P(0) > + Ons (5-30) 
P, | h*0*P(0) 
Ps h505P(0) 


Recall that the matrix in equation (5-30) is denoted [AA] as before. 


After multiplying each side of the equation by caay"!, the result becomes, 


(0) Py 


h2D2P(0) Pp 
h303P(0) >= (Ma) Ps) = One (5-31) 
h40*P(0) | P, 
hSD5P(0) Ps 


Now P(0) is approximated by the coefficients in the first row of caay 
and D2P(0) is approximated by the coefficients in the second row of 
caay > divided by h2, both expressed in terms of P, through P,. Notice 
that the approximation for D2P(0) has truncation error Oh* and P(0) has 
truncation error Oh®, which is acceptable in this analysis. 

The approximation for D2Q(0) is similar except that Q(0) will appear 
explicitly in the approximations with DQ(0) = 0. Equation (5-10) now 


becomes, 
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{Q(0)} + 


/2029(0) . 
h?D7Q(0) | 
h*0*Q(0) 

[AA] + Oh® 
6x6 h805Q(0) 
h*p*Q(0) 
h707Q(0) 

(5-32) 


After rearranging, equation (5-32) is multiplied through by caay and 


becomes, 


h202Q(0) 
h303Q(0) 
h4D*Q(0) 
h§D5Q(0) 
h®0%Q(0) 
h707Q(0) 


= aay} 
6x6 


Q:- Q(0) 
Qe- Q(0) 
Qs- Q(0) 
Q- Q(0) 
Qs- Q(0) 
Qe- Q(0) 


- Oh & (5°33) 


The coefficients for 02Q(0) are taken from the first row of caay + and 


divided by h?. 


02Q(0) is now expressed in terms of Q(0) and Q, through 


Qg and will have a truncation error Oh®. 


The two additional boundary condition equations involving y can now 


be solved in the eigensystem of equations with H(0) and Q(0) appearing 


explicitly as the two required unknowns. It is convenient to express 


the [Cg] and [Dg] matrices of equation (5-29) in slightly different 


form, 
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[Cy] I HCO) as [Cy] {H(0)} 
and | (5-34) 


[09] at HCO) as {DeT {H(0)} 


where [C,*] and (0.*] become 2 X 1 column matrices, 


(= - 4+ zat) “4 
(C,*] = and [D9*] = (5-35) 
~202 + ine 1 


Substituting equations (5-35) into equation (5-29) yields, 


36. (07?) P(0) 
- Fe [Cr] + 2ia[Cg] + ia{Cg*] {H(0)} 
02Q(0) Q(0) 
(5-36) 


P(0) 
= y {2 [0g] + a2 (D9*] {H(0)} 
Q(0) 


The solution of the eigensystem of 2N + 2 equations in 2N + 2 un- 
knowns is now accomplished in the form of equation (5-3). The appear- 
ance of the [A] and [B] matrices are identical to those for n 2 6 in 
Figures 5-5 and 5-6 except for the addition of two columns to accommo- 
date the parameters H(0) and Q(0), and two rows to accommodate the pair 
of boundary equations in (5-36). Nonzero elements are present in the 


{A] and [8] matrices for n = 1 as shown in Figures 5-7 and 5-8. 
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Figure 5-7. [A] Matrix for n=1 
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Figure 5-8. [(B] Matrix for n=1 
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Derivation of the noncentral finite difference coefficients for the 
derivatives of P and Q at the wall are straightforward as shown for the 
other angular wave numbers. The columns of the [AA] and [CC] matrices 
are eliminated based on the fact that parameters P(1), Q(1), and 0Q(1) 
are all expressed in terms of H(0), as seen in equations (3-18c). By 
imposing the boundary conditions for Q at the wall, 


Q(1) = - H(0) and = DQ(1) = 2 HCO) 


equation (5-17) becomes, 


Qi + (1-h) HCO) | 
Qo + (1-3h) H(0) 
Qs + (1-5h) H(0) 
Qq + (1-7h) HCO) 
Qs + (1-9h) HCO) 
Qe + (1-11h) H(0) 


{H(0)} + (CC) (AAJ 2 
4x6 6x6 


(5-38 


The coefficients in the resultant 4 X 6 matrix of equation (5-38) 
are in fact the same coefficients computed for n 2 6 at the wall, and 
recall that they appear in reverse order when used at the wall. The 
details of the complete derivation of the noncentral finite difference 


coefficients near the boundaries jis lengthy and not included here. 
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D. COMPUTER PROGRAM USEAGE 


The solutions for the discrete transformed vorticity transport 
equations for the three angular wave numbers investigated are ac- 
Be oaks complished in the main investigative computer programs. Three separate 
programs were written rather than developing one program for the general 
case. The rationale for this is as follows: the changes of variables 
for each angular wave number n $ 6 are different, the boundary condi- 
tions are different, the equations uncouple for n = 0, additional para- 
meters are introduced for some, and n = 1 has a coupled boundary 
equation that contains y that must be dealt with separately in the setup 
of the problem. It would be practical to write a general program for 
angular wave numbers n 2 6, however. 

For exch specific angular wave number analyzed, the main computer 
program is divided into eight parts. Their functions are summarized 
below: 


Part I The central finite difference coefficient for the derivatives 
of P and Q are computed. 


Part II The noncentral finite difference coefficients for the deriva- 
tives of P and Q at the axis and wall] are read in as data. 
Problem variables, Re and a, are read in as data. 


Part III Coefficients of the vorticity transport matrix equations are 
computed along the radial mesh. 


om Part IV The [A] and (B] matrix elements are computed. 

; Part V The eigenvalues and respective eigenvectors are computed in 
mA sy the subroutine EIGZC and the least stable eigenvalue is 
A determined. 


: Part VI Verification that the least stable eigenvalue and correspond- 
1 ing eigenvector satisfy the governing equations. 
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Part VII The normalized axial perturbation velocity is computed from 
the least stable eigenvalue and its eigenvector. 


Part VIII The normalized axial perturbation velocity is plotted versus 
normalized pipe radius utilizing the Versatec plotter. 

The numerical analysis of the pipe flow stability problem jis cen- 
tered around the subroutine EIGZC once the problem is properly set up. 
The eigensystem of equation (5-3) is solved by this subroutine which 
returns a set of N, 2N + 2, or 2N eigenvalues for n = 0, 1, or 6, 
respectively, and a corresponding set of normalized eigenvectors. Each 
eigenvalue and its eigenvector represent a solution of the governing 
equations. The primary objective of the main program and of this in- 
vestigation is to determine the least stable eigenvalue which indicates 
flow characteristics. Recall that if ¥p is positive, the flow is 
unstable. For each value of a and Re chosen, the maximum algebraic 
value of yp is found, denoted as y*, to determine if instabilities exist 
in the flow field. The corresponding eigenvector of y*, denoted as 
{X*}, is used to compute the incremental axial perturbation velocity 
along the radius. Plots of the normalized perturbation velocity versus 
normalized pipe radius have no significance in determining stability or 
instability of the flow but are indicators of the general nature and 
behavior of the flow field perturbations. Another independent program 
was written to plot Reynolds number stability contours using data taken 
from program results for n = 1 and 6. The subroutines used to make 
plots on Versatec plotter are PLOTG, standard Versatec plotting func- 


tions and TITLE1 to print titles and legends on the pilots. 
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The numerical accuracy of the solution was determined by substi- 
tuting y* and its eigenvector into the governing equations. For this 
purpose, the residual error vector {ce} was computed by equation (5-39) 


and the magnitude of the error was examined. 


fe} = [A] {X*} - y* (BJ {X*} (5-39) 


The Fortran subroutine EIGZC jis contained in the International 
Mathematical and Statistical Library (IMSL) available in most large 
computer system libraries utilizing Fortran compilers. Specifically, 
the numerical analysis of the vorticity transport equations in the form 
of equation (5-3) was computed on the IBM 370/3033 (assigned processor) 
system with the Fortran H Extended compiler at the U.S. Naval Post- 
graduate School. All calculations were performed in full double pre- 
cision format to improve the accuracy of the solution by minimizing 
truncation and round off errors inherently introduced by digital compu- 


tations. The main investigative programs and other auxiliary programs 


were all written in Fortran programming language to be compatible with 


all versions o. Fortran compilers. The computer program listings with 


applicable data sets are included after the appendices in this paper. 
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VI. RESULTS 
A. RESULTS FOR n = 0 


Results of the pipe flow stability problem are obtained from the 
solution of the linearized vorticity transport equations by the methods 
previously described. A generalized solution of the governing differen- 
tial equations for fully developed, three-dimensional pipe flow has been 
accomplished for three angular wave numbers, n = 0, 1, and 6. ‘The 
characteristics of the flow field perturbations for each angular wave 
number investigated are established by the parameters Re and a. The 
purpose of this numerical investigation was to determine flow stability 
at various Reynolds numbers and axial wave numbers for each angular wave 
number. Recall that the perturbation velocity vector Vv is defined in 
equation (2-27) as the curl of the velocity vector potential W. This 
function is expressed in terms of G(r) and H(r) and a complex exponen- 
tial form e* in equations (2-28) and (2-29). After appropriate changes 
of variables were introduced, solutions were determined from the eigen- 
system of equations which consists of a set of eigenvalues and their 
corresponding normalized eigenvectors for fixed values of Re anda. The 
resulting perturbation quantities now appear in terms of the eigenfunc- 
tions P(r) and Q(r) and are represented in discrete form in the eigen- 
vector {X} as P; and Q;; where i = 1, 2, 3..... N. As can be seen from 
equation (2-30), the real part of the complex eigenvalue Yp will deter- 


mine the growth or decay rate in time of the perturbation. Positive 
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r ‘values of Yp indicate that the flow is unstable. The eigenvalue whose 

( real part has the largest algebraic value, denoted as y*, is the least 
; stable and is reported in these results to show stability trends. 

Gawain (Ref. 13] reported the numerical results for angular wave 

number n = 0 for the solution of the governing equation in his paper. 

Those numerical results which incorporated a purely imaginary axial wave 

number and the advancements in the derivation of the boundary conditions 

" at the pipe axis were duplicated in this investigation. They are 

essentially the same results except for variations in the fourth decimal 

place of the value of y*. This is due to the improved accuracy of the 

i : solution as a result of the numerical methods used here. The primary 


emphasis of this investigation is centered on determining the onset of 


flow instabilities at a theoretical critical Reynolds number near 1150. 
( However, solutions were obtained at other Reynolds numbers to show 

stability trends and are presented in tabular and plotted form. 
: Stability data for angular wave number n = 0 for various combin- 
ations of Re and a are given in Table 6-1. The least stable eigenvalue 
solution of the governing equation jis based on a computational mesh of 
N= 50. Table entries which contain dashes indicate that the solution 
of the governing equations was not sufficiently accurate to be included 
as valid data. Two facts become very clear in examining Table 6-1. 
First of all, the flow is stable at all Re anda for n=0. Secondly, 
there is an asymptotic trend towards neutral stability for increasing 
Reynolds numbers. This second result agrees with experimental observa- 


tions in that increasing Reynolds number has a destabilizing effect on 


the flow. It can also be seen that increasing the axial wave number has 
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a stabilizing effect on the flow. Normalized axial perturbation velo- 
city plots versus pipe radius were made at various Re and a to examine 
the behavior of the perturbation quantities corresponding to the least 
stable eigenvalue solution. Figures 6-1 through 6-9, on the following 
pages, are representative of the behavior of the axial perturbation 
velocity. The activity is concentrated near the pipe axis for n = 0 and 
dies out quickly as r approaches 1.0. Additionally, the perturbation 
velocity plots were made in order to determine if mesh fineness was 
sufficient to obtain an accurate solution of the equations. Trial 
calculations were made for uniform mesh sizes up to 100 where the 


numerical results did not vary appreciably from those reported here. 
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B. RESULTS FOR n = 1 


The numerical results for angular wave number n = 1 indicate that 
flow instabilities indeed exist. There are four factors that may 
possibly account for the flow instabilities; (1) the complicated nature 
of the coupled governing equations for n = 1, (2) the rigorous way in 
which the boundary conditions at the axis were derived and enforced, 
(3) the introduction of additional parameters into the eigensystem, and 
(4) the appearance of the eigenvalue in a pair of the boundary condition 
equations. Table 6-2 summarizes the stability data for n = 1 over a 
wide range of Reynolds numbers and axial wave numbers. Once again the 
dashed table entries indicate inaccurate results. Notice that a mesh 
size of n = 48 was used for this case. This was necessary because of 
programming constraints that required the eigensystem to be less than or 
equal to 100 X 100. It can be seen that flow instabilities exist at all 
Reynolds numbers. This occurs coincidentally with smal] values of a. 
More detailed stability data is presented in Table 6-3 for axial wave 


numbers 0 S$ as 1. It is clear here that the flow becomes more unstable 


for very small values of a <1. The stability data presented in Tables 


6-2 and 6-3 is more easily interpreted as Reynolds number contours in 


be 
fis 


the plots of GAMMA*® vs. alpha. Figure 6-10 shows Reynolds number con- 


tours over a large range of alpha. It can be clearly seen that the flow 


Woo 


is unstable at all Reynolds numbers for small values of axial wave 
f numbers. Figure 6-11 shows a more detailed plot of the Reynolds number 
: contours over a small range of a. Here it can be seen that the flow 
f becomes more unstable for decreasing values of Reynolds numbers! It can 
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also be seen here that the flow becomes even more unstable for small 
values of a only slightly larger than zero. The reason for these 
results cannot be readily explained but it is obvious that they are in 
direct conflict with observed experimental results. When a > 1, however, 
the contours show an asymptotic trend toward neutral stability for 
increasing Reynolds numbers; in other words, increasing Reynolds number 
has a destabilizing effect on the flow which agrees with experimental 
results. A slightly different method of presenting the data can be seen 
in Figure 6-12, as alpha contours in the plot of GAMMA* vs. Reynolds 
number. This plot further emphasizes the results discussed above; the 
flow is unstable at all Reynolds numbers for smal] alpha, increasing 
Reynolds number has a destabilizing effect on the flow and increasing 
axial wave number has a stabilizing effect. 

Normalized perturbation velocity plots were also produced for this 
case and show some interesting trends. The behavior of the axial per- 
turbation velocity for selected values of Re and a appear in Figures 
6-13 through 6-23 on the following pages. Upon close examination of the 
perturbation velocity plots, it becomes apparent that the computational 
mesh is not sufficiently fine to adequately approximate the velocity 
function as a is increased. This problem will be discussed shortly. 

It is interesting to note that the perturbation activity shifts from 
the pipe wall to the axis for increasing values of the axial wave number 
while Reynolds number remains constant. At small values of a, the 
activity originates at the wall. At larger values of a, the activity 
shifts to the axis, for the most part. At even larger values of a, the 


activity moves away from the axis to the region between the axis and 
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wall. However, some degree of activity remains at the wall for the last 


two cases. The change in the perturbation activity as described above 
can be seen in Figures 6-14 through 6-16 for Re = 575 and in Figures 
6-19 through 6-21 for Re = 2300. The shift of perturbation activity to 
the axis with increasing axial wave number deserves comment here. 
Schlichting [Ref. 16] observed that the transition to turbulence is 
“characterized by the ........ appearance of self-sustaining turbulent 
flashes with emanate from fluid layers near the wal) along the tube." 
Therefore, it appears that unstable flows are associated with pertur- 
bation activity near the wal]. As the activity shifts away from the 
wall with increasing axial wave number, the flow becomes stable. Addi- 
tionally, after the activity moves to the axis, some degree of activity 
remains at the wall. While the solution appears to remain stable, it is 


recommended that this peculiarity be examined further. The use of a 


finer mesh to resolve this activity at the wall would be useful here. 
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C. RESULTS FOR n = 6 


Results for angular wave number n = 6 are similar to those for n = 0 
in that the flow is stable at all Reynolds numbers. Stability data is 
presented in Tables 6-4 and 6-5 where it can be determined that this is 
the case. The data is once again presented as Reynolds number contours 
in the plots of GAMMA* vs. alpha in Figures 6-24 and 6-25 and as alpha 
contours in the plot of GAMMA* vs. Reynolds number in Figure 6-26. The 
trends are clear; the flow is stable at all Reynolds numbers, increasing 
Reynolds number has a destabilizing effect on the flow and increasing 
the axial wave number has a stabilizing effect. 

Figures 6-27 through 6-36 are representative perturbation velocity 
plots for n = 6 at several Reynolds numbers and axial wave numbers. It 
can be seen in Figures 6-31 and 6-32 that the perturbation activity 


shifts from the wall to the axis for increasing a. 
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Ee D. NUMERICAL ACCURACY 
A. Numerical accuracy of the solutions of the governing vorticity 
transport equations was verified by substituting the least stable eigen- 
value and its corresponding eigenvector into the eigensystem relation 
given in equation (5-39). The error vector was examined to determine 
the residual error of the solution. For the most part the error was of 


the order 10°22 to 1072 


» a most satisfactory result. 

Upon closer examination of the error vector and the stability trends 
in the plots, it became evident that the accuracy of the numerical 
solution for the governing equations was not only dependent on the 
Reynolds number as previously mentioned but on the axial wave number as 
well. Plots of the perturbation velocity vs. radius for very large 
axial wave numbers are contained in Figures 6-37 through 6-39. The 
three figures, one for each angular wave number investigated, are 
examples of inaccurate solutions. Dependence of the solution on Rey- 


nolds number and axial wave number, became more evident as Re and a were 


increased. Rapid excursions of the perturbation velocity over a smal] 


wo 
fn 
na 


portion of the pipe radius indicate that the computational mesh was not 


sufficiently fine to approximate the perturbation velocity function in 


q aoe 


the region of activity. This supposition can also be supported by 
examining the error vector where the residual error was found to be of 
the order 10°! or worse, in the region of perturbation activity. Recal] 
for n = 1 that some of the activity remained at the wall as a was in- 


creased. Refer to Figures 6-20 through 6-23. The residual error was 


wry wpe eee peop 
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nominally of the order 10°? in this region. Therefore, the true repre- 
sentation of the perturbation velocity is directly dependent upon an 
accurate numerical solution of the governing equations. These inac- 


curacies can be resolved by using a nonuniform mesh in the region of 


interest as suggested by Arnold. It is suspected, however, that after a 


thorough analysis of the effects of the problem parameters and mesh size 
on the solution have been investigated, only small values of the axial 
wave number below about 10.0 and Reynolds numbers in the vicinity of 
1150 will be of primary interest. It is felt by this investigator that 
a mesh size of N = 50 and the wide range of problem parameters used here 
yielded solutions that satisfied the governing vorticity transport 


equations to a high degree and were sufficient to show stability trends. 
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VII. 


RECOMMENDATIONS AND CONCLUSIONS 


The problem of predicting the transition from laminar to turbulent 
flow of incompressible fluids in circular pipes has been pursued for 
nearly 100 years and has been presented here in a straightforward manner. 
The governing equations for the dynamics of this problem consist of the 
Navier-Stokes equation and the continuity equation. After taking the 
curl of the Navier-Stokes equation and applying the continuity principle, 
the linearized vorticity transport vector equation is obtained. The 
velocity vector potential function, which contains functions of r and a 
complex exponential factor, was then introduced. The resulting two 
linearily independent equations form the basis for the solution of the 
problem in the form of the vorticity transport matrix equation. 

A valid solution of the governing equations was made possible be- 
cause of two factors; (1) the advancements in the linearized theory by 
introducing a purely imaginary exponential form of the axial wave number, 


and (2) the rigorous method in which the boundary conditions at the pipe 


axis were derived and those constraints enforced. After an appropriate 
change of variables was performed, the governing equations were trans- 
formed into an equivalent pair of coupled, homogeneous differential 
equations that were solved in eigenvalue problem format. Accurate 
solutions of the governing equations were achieved by the refinement and 
combination of several standard methods of numerical analysis. These 


methods include; (1) the approximation of the governing equations along 


a computational radial mesh developed by the half-station method, (2) 


the use of central and noncentral finite difference approximations 
having consistent fourth order truncation error, and (3) the use of 
complex, double precision number representations. 


It was this author's intention from the outset to make this work as 


_ complete as possible. This was accomplished by including the develop- 


ment of the theory, the derivation of the boundary conditions, and a 
detailed explaination of the numerical methods used as well as obtaining 
solutions to the governing equations. This will enable future investi- 
gators to use this thesis as a starting point and as a single source 
reference for follow-on studies of this subject. 

The conclusions of this numerical analysis follow directly from 
the results of the pipe flow stability problem which were discussed 
thoroughly in the previous chapter. 

The most significant findings are: 

(1) for n = 0, the flow is stable at all Re and a. 

(2) for n = 1, the flow appears to be unstable at all Re for small 

values of a. This result is unfortunately in direct conflict 


with experimental observations. 


(3) for n = 1, unstable eigenvalues are associated with perturbation 
activity at the wall. 


(4) for n = 6, the flow is stable at all Re and a. 
The following observations are also of interest: 


(5) for n = 0, 1, and 6, increasing Reynolds number has a destabi- 
lizing effect on the flow. 


(6) for n = 0, 1, and 6, increasing axial wave number has a stabi- 
lizing effect on the flow for fixed Re. 


(7) for n = 1 and 6, increasing axial wave number causes a shift in 
the perturbation activity from the wall to the axis for fixed Re. 


The results for angular wave numbers n = 1 and 6 are believed to be new. 
Conclusions (1), (5), (6), and (7) are in agreement with experimental 
results or previous theoretical investigations. 

It is this author's opinion that the numerical analysis of the 
vorticity transport equations is now complete for angular wave numbers 
n=0 and 6. It is felt, however, that the numerical results for 
angular wave number n = 1 warrant further study. The unresolved paradox 
is quite puzzeling in that for n = 1 the flow appears to be unstable at 
all Reynolds numbers. Continued numerical investigations for n = 1 and 
follow-on studies of other angular wave numbers should prove fruitful in 
establishing a theory which predicts the transition from laminar to 


turbulent flow of fluids in pipes. 
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APPENDIX A 
COEFFICIENTS OF THE VORTICITY TRANSPORT EQUATIONS 
The generalized coefficient matrices which appear in the basic 
vorticity transport equations (2-40) are defined below. 
( 0 0 
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APPENDIX B 


COEFFICIENTS OF THE TRANSFORMED VORTICITY TRANSPORT 
EQUATIONS FOR n = 0, 1, and 6 


A. COEFFICIENTS FOR n = 0 


Since the vorticity transport equations uncouple for n = 0, only the 
solution for the eigenfunction Q(r) is sought. Therefore the details of 
the complete coefficient matrices are not required. The coefficients 
which appear here represent the matrix element (2,2) and correspond to 


the primed quantities of equation (2-43). 


Mo =- a5 (B-1) 
Ms = - £- (B-2) 
Ms = a5 ( 34 202+ + 2ia (r - r3) (B-3) 
Mi = (3 + ea?) + Gia (1 - r2) (B-4) 
Mo = - ar + 2ia? (r3 - vr) (B-5) 
No =-r (8-6) 
Ni; = - 3 (B-7) 
No = a?r (B-8) 
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These coefficients also appear and are computed in part III of the main 


ie investigative computer program for n = 0. 
B. COEFFICIENTS for n= 1 


Since the vorticity transport equations do not uncouple for n = 1, 

the eigenfunctions P(r) and Q{r) must be solved for simultaneously. 

With the introduction of the additional parameter as a result of the 
Nf change of variables, H(0) appears explicitly as an unknown in the system 
of equations. This requires the addition of two coefficient matrices to ! 

equation (2-40). The 2 X 1 column matrix [Ms] appears on the left side 

a of equation (2-47) and one 2 X 1 column matrix [Nz] appears on the right 

side of equation (2-47). The complete coefficient matrices for n= 1 


are defined below. 


0 0 
[Ma] = 5 | (8-9) 
0 7 Re 
a ae 
[M3] = Re (B-10) 
oe Woe 
Re 
e (Mz) = (B-11) 
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a2 
These coefficients also appear and are computed in Part III of the main 
investigative computer program for n = 1. 
C. COEFFICIENTS for n = 6 
The vorticity transport equations for n = 6 remain coupled. There 
q are no additional parameters introduced with the change of variables to 
P(r) and Q(r). The system of equations take the form of equation (2-40) 
once again. The complete coefficient matrices for n = 6 are defined 
a below. 
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These coefficients also appear and are computed in part III of the main 


investigative computer program for n = 6. 
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APPENDIX C 
COEFFICIENTS OF THE BOUNDARY EQUATIONS AT THE AXIS 


The matrices which appear in the boundary equations at the axis, 


equations (3-4) through (3-8) are defined below. 


- n2 - in 
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+4 in - (n2 + 3) 
= - n2 - 2 in 
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% +4in  - (n2 + 4) 
7 ; | 
: Bodh 
. [C3] = (C-3) 
' +2in  - § (n? + 3) 
: 
M . 
? - }2n2 - (2n2 + 1) fe, - in (2 - za) 
: [C4] = (C-4) 
- +4 in (1 - 32) -2(n2 +1) (1 - 2] 
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fee + 202 +2 = 202] 


Column matrices used in equation (5-36) for n = 1. 


[Co*] = 


[Dy*] = 
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APPENDIX D 
SPECIAL CONDITIONS AT THE AXIS 


The determinates are formed from the corresponding matrices as 
defined in Appendix °. For angular wave numbers n 2 6, all of the 


determinates become non-singular. 
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